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Abstract 

^ The success of an infectious disease to invade a population is strongly 

'"5 controlled by the population's specific connectivity structure. Here a net- 

^\ work model is presented as an aid in understanding the role of social be- 

havior and heterogeneous connectivity in determining the spatio-temporal 
patterns of disease dynamics. We explore the controversial origins of long- 
term recurrent oscillations believed to be characteristic to diseases that 
have a period of temporary immunity after infection. In particular, we 
focus on sexually transmitted diseases such as syphilis, where this contro- 
versy is currently under review. Although temporary immunity plays a 
key role, it is found that in realistic small-world networks, the social and 
sexual behavior of individuals also has great influence in generating long- 
term cycles. The model generates circular waves of infection with unusual 
spatial dynamics that depend on focal areas that act as pacemakers in 
the population. Eradication of the disease can be efficiently achieved by 
eliminating the pacemakers with a targeted vaccination scheme. A simple 
difference equation model is derived, that captures the infection dynam- 
ics of the network model and gives insights into their origins and their 
• eradication through vaccination. 

1^ 

Developing strategies for controlling the dynamics of epidemics as they spread 
r through complex population networks is now an issue of great concern HI HI |3J 

m [3 El [3 [Hi ini 113 • Future progress depends on gaining a better theoretical 
^ understanding of the spatial dynamics of disease spread, including the effects of 

a population's social contact structure and its network topology [51 [51 HI 151 lllj. 
%N Here we show how these factors control epidemic spread and, in the process, 

formulate a novel aggregated targeted vaccination scheme. 

We are particularly interested in diseases that confer temporary immunity 
to individuals after recovery from infection. This is typical for diseases such as 
pertussis, influenza, hRSV and some sexually transmitted diseases (STD's) as 
syphilis. In terms of population dynamics, the temporary immunity is under- 
stood to give rise to recurrent epidemic oscillations [12] that can have a period 
of several years for pertussis [7] and certain strains of influenza [T3] , to decadal 
oscillations in the case of syphilis [TUITn]. In simple terms, the epidemic cycles 
arise due to a delayed " SIRS" process in which Susceptible individuals become 
Infected, Recover with temporary immunity, but then eventually return to the 
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Susceptible pool after a time-delay when immunity wears off. The loss of im- 
munity allows the susceptibles in the population to gradually build up until 
sufficient in number to fuel the next disease outbreak. 

Crassly et al. [Hj suggested that the oscillations seen in long-term syphilis 
data-sets from the US stem from the temporary immunity of this disease. Their 
argument is buffered by the fact that gonorrhea, which lacks temporary im- 
munity, fails to show the same strong cycles in long-term datasets. This view, 
however, is controversial and the CDC [TB] has countered that trends in US 
syphilis epidemiology follow parallel changes in population-wide high-risk sexual 
behavior (see also [T5] ) . Most likely it is the combined presence of temporary im- 
munity and social behavior that is responsible for the recurrent waves of syphilis 
epidemics. The modeling approach described here allows us to investigate and 
assess the impact of these different but important factors. 

Complex networks (or graphs) provide an important means for investigating 
the effects of social behavior in population models of disease spread. Individu- 
als are represented as nodes of a graph and edges are placed between any two 
individuals should there be an infection route between them [SJIUI^. A random 
Erdos Renyi network is formed if there is an equal probability g of a connection 
between any two individuals |17j . A regular and tightly clustered network struc- 
ture is obtained if an individual is only able to infect his/her nearest neighbors. 
The random and clustered-regular graphs might be considered as two cndpoints 
of a spectrum. Watts and Strogatz [TH] developed a scheme that allows con- 
struction of networks that interpolate anywhere between these two endpoints. 
This is achieved by introducing a proportion of p random "short-cuts" between 
nodes in a regular graph. Only relatively few short cuts are required {p < 0.1) 
to create "small world" networks that have the often realistic qualities of both 
a high degree of clustering, and at the same time relatively high overall network 
connectivity introduced by long-range connections (i.e. via short-cuts). 

When considering the population dynamics of STD's it is important to take 
into account that some individuals spread the disease to a much greater extent 
than others. In this way, social behavior and sexual promiscuity governs the 
heterogeneity of the contact structure in the population. This contrasts with 
standard mean field differential equation models which are based on the as- 
sumption of a "randomly mixing" population and lack a heterogeneous contact 
structure. However, there is no unanimous agreement on how the contact struc- 
ture of the network should be fixed. Barabasi et al. [20l [21] have argued that 
"scale free" networks, whose nodes have a power law connectivity distribution, 
are the most appropriate for STD's. Lloyd and May |TI], on the other hand, 
suggest that such a formulation is unnecessarily exaggerated. We follow Eames 
and Keeling [3J 0] who use a small world model as a first approximation for 
STD's. Their model assumes that STD's are generally transmitted locally but 
long-range infection pathways exist and are important in spreading the disease 
through the population network. Moreover, the small world approach creates 
heterogeneity in the connectivity distribution with most individuals having sev- 
eral connections, but some being more connected than others. 

The small-world formulation is used here to study the spatio-temporal SIRS 
dynamics of recurrent diseases with temporal immunity. We first describe the 
network model and its spatio-temporal dynamics. For representative parame- 
ters the model exhibits expanding circular waves of infection, some of which are 
generated by unusual "pacemaker centers". These we study in detail; the im- 
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Figure 1: Spatial SIRS model simulation illustrating disease dynamics and ag- 
gregated vaccination. Parameters (as in [III): k — 3,tj — A,tr — 9,p — 
0.02, q — 0.2. (a) t — 10. Circular waves of infected individuals (red) spread 
through a population of susceptible (white) and recovered/immune (green) in- 
dividuals located in a 200 x 200 lattice, (b) t — 40. Data-analysis has identified 
two periodically reappearing pacemaker centers (red rings) of infected individ- 
uals, (c) t = 90. Aggregated vaccination of all individuals (blue) located in 
proximity to the pacemaker centers - comprising only 18% of the entire pop- 
ulation, (d) t = 100. Disease rapidly brought to extinction (all red infectives 
eliminated) in absence of other pacemaker centres. 



portant role of such "pacemakers" suggests a practical disease control strategy 
based on targeted vaccination. We show that by vaccinating or quarantining the 
regions around pacemaker centers, the disease can be eradicated. The vaccina- 
tion scheme is tested on various more realistic modifications of the basic model. 
We then formulate a very simple difi^erence equation model that captures some 
of the main features of the heterogeneous network. 

1 The network SIRS model 

The network model is based on a 2— dimensional lattice of individuals (/nodes) 
whose connectivity p can be preassigned. Each node on the lattice is occu- 
pied and is connected to k nearest neighbors oriented in each of four directions 
(North, South, East, West, with diagonal connections excluded). That is, each 
node is initially connected to K = Ak nearest neighbors, with fc = 3 unless oth- 
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erwise specified. The horizontal and the vertical edges of the lattice are "glued" 
together creating a 2— torus. Then, with a probability p, each of the K nearest 
neighbors of each of the edges in the lattice is randomly rewired to an arbitrary 
node. These rewired connections, or "short-cuts" [HI US], may extend to far 
regions of the network. 

The parameter p controls the population's connectivity structure: p — Q 
corresponds to nearest neighbor contacts only, and where clustering is at its 
maximum; small p in the range < p < 0.1 corresponds to a "small world" 
network (each individual has a certain amount of nearest neighbor contacts + a 
small proportion of distant contacts, "short-cuts"); large p > 0.4 is qualitatively 
equivalent to a randomly mixing population with minimal clustering |19j . As 
p is the probability of a short-cut, it may be viewed as an index of population 
mobility. Alternatively, it may be interpreted as an index of social behavior 
such as sexual promiscuity in the case of STD's, given the manner in which it 
controls overall network connectivity and clustering |18j . 

Disease dynamics follow the classical Susceptible-Infectious-Recovered-Susceptible 
(SIRS) formulation [H UHl H E] with Susceptible individuals (S") having a 
probability q of becoming infected when linked to an Infected individual (/). 
Infected individuals eventually recover from the disease after a fixed time pe- 
riod, Tj, and are conferred temporary immunity. After a time period of th time 
units, immunity wears off and Recovered individuals (R) return once again to 
the Susceptible pool (S) closing the SIRS loop. 

This is implemented on the network using a 2— dimensional cellular automata 
(CA), SIRS spatial model. At time t, an individual at the (j,j)'th location of 
the lattice has the state Xijit) which is cither I or R. The model is based on 
the following transition rules: 



Xi^j{t) el — > x^^j{t + 1) e / ^ • • • — > Xij{t + Tj) e R 
Xij{t) e R^ Xi^j{t + 1) e i? ^ • • • — > Xi^j{t + tb) e S. 

Infections are transmitted to susceptible individuals with a probability q, if 
they are connected to an infective via a nearest neighbor or a short-cut. Thus 
the probabihty that a susceptible becomes infected is 1 — (1 — q)'''"' , where /cjnf 
is the total number of infected neighbors of the individual, be they connected 
via nearest neighbors or via short-cuts. The proportion of S(t),I{t) and R{t) 
individuals are calculated over the lattice and their dynamics are followed as a 
function of time. 

To help fix ideas, we focus on two representative parameter settings: (a) 
parameters used in a general theoretical model taken from [19]. The infectious 
period is fixed at r/ = 4 time units and a recovery period of r/j = 9 time units 
(see simulations in Fig.[T]); (b) parameters associated with syphilis epidemics as 
based closely on the study of Crassly et al. [T3] (see simulations in Fig. |2| . In 
the latter case r/ = 1 time unit, which is taken to correspond to half a year; the 
recovery time varies randomly and uniformly in the range S {8, 9, 10, 16} 
time units corresponding to a period of 4 — 8 years of immunity. To add realism, 
several other features are also incorporated. A birth-death process is introduced 
at rate /i, indicating the proportion of births/death per time step with appro- 
priate network rewiring. In Fig. [2| = 0.01 per half-year time step, which is 
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prob. 1 — (1 — g) 
otherwise 
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Figure 2: Model simulation of syphilis dynamics white - susceptibles, red - 
infected, green - recovered, k — 3, tj — 6 months, tt^, = 4 — 8 years, p — 0.065, 
q — 0.36, ip = 0.2 p.y., fi — 0.02 p.y. (a) Three pacemakers are seen at time 
t — 12.5 years, (b) Pacemaker areas are vaccinated with random spread of 
86% of each pacemaker area on average, at t = 41 (20.5 years), resulting in 
vaccinating 14.5% of the population, (c) Disease eradication at t = 58 (28y.). 
(d) Time series of the proportion of infectives per 0.5y. The period of oscillation 
is r w lOy. The vaccination time is indicated by 'o', after which the disease 
undergoes another smaller peak, and reaches complete extinction within less 
than one period of the disease (8.5 years). 



equivalent to a rate of fi = 0.02 per year. Provision is made for the possibility 
that a proportion (p of individuals fail to gain immunity after infection (similar 
to [H]). Thus If = corresponds to all nodes passing through an SIRS loop 
while (fi = 1 corresponds to all nodes exhibiting SIS dynamics (no individual can 
acquire immunity). In Fig. |2]a proportion of = 0.2 of the population (20%) 
gaining no immunity per year. 

1.1 Recurrent circular waves and pacemaker centers 

We have found that in the small world regime, models of type ([T]) exhibit 
concentric waves which give rise to unusual "pacemaker centers." Indeed, for 
0.001 < p < 0.15, the model (jlj exhibits spatial oscillations with expanding cir- 
cular waves of infection traveling through the lattice (see Fig.jTJi). Some of these 
waves are recurrent, both spatially and temporally. The latter are generated 
by pacemakers (Figs. Qp and ^a) that form at connectivity centers - localized 
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areas denser in short-cuts. The waves grow in size about the pacemakers as 
the infection spreads radially. When infected individuals recover, the interior 
of the growing wave boundary becomes a fresh pool of susceptible individuals. 
At the end of the cycle, a distant infectee "short-cuts" through the network to 
reinfect the wave's focal pacemaker, enabling it to perpetuate. The cycle allows 
recurrent spatial waves to propagate with a fixed period, T. For the syphilis 
parameters (see Fig.[2]legend), T w 10 years, as observed in US syphilis datasets 
[Til [53] . Similar spatial wave patterns had been observed in a number of bio- 
logical contexts including epidemiology [24J, ecology [25], neural networks [26] 
and theoretical studies of excitable systems [?7] . 

These pacemakers follow a pattern formation mechanism. We observe that 
there is a minimal necessary amount of short-cuts needed for the creation of a 
pacemaker center. In addition, within the small-world range, a pacemaker wave 
region has more short-cuts than a non-pacemaker wave region. 

The initiation of a pacemaker also requires that the infection is able to spread 
both in the horizontal and in the vertical directions. That is, the development of 
the infection from an initial state should progress in the two directions spanning 
a plane (i.e., in an X or an L shape). This condition follows the rule that 
heterogeneities are needed for creation of spirals in excitable media (see [27]). 
For these reasons, having an aggregation of short-cuts in a small region enhances 
the likelihood of creating a pacemaker. On the other hand, having too large 
an aggregation of short-cuts in a localized area of the lattice results in the 
opposite effect - disease extinction will occur in this localized area due to a 
synchronization effect (see [2S] and below). 

1.2 A Targeted Vaccination Scheme 

The important role of "pacemaker centers" suggests a practical control strategy. 
We have found that by vaccinating or quarantining the regions surrounding 
pacemakers, the disease can be usually brought to a complete extinction within 
one period of the disease (in some cases, depending on the refinement of the 
vaccination algorithm and specific parameter values, two vaccination pulses are 
required). Thus, rather than the conventional scheme of immunizing some 85% 
of the population to achieve herd immunity [1] , it is only necessary to vaccinate 
groups enclosing the pacemakers. This requires vaccination of some 10%— 30% of 
the population (depending on the specific application and algorithm refinement; 
~ 20% in most cases). Fig. [ijand Fig. |2]show spatial snapshots of an infected 
population upon application of the vaccination scheme. In the first frames 
(subfigs. la,b and 2a) the characteristic circular waves of infection (red) are 
seen. Vaccination around the pacemakers leads to complete eradication of the 
disease. Pacemakers have such large impact on the spatial dynamics that they 
are relatively easy to detect using a simple threshold algorithm that identifies 
recurring aggregations of infected individuals. Once a pacemaker is identified, a 
small region enclosing the pacemaker is marked out and vaccinated by effectively 
removing these nodes from the simulation (the blue rectangles in figures [lj3,d 
and|2|3,c). 

For the theoretical values of [19^ used in the simulations of Fig. [T| the ratio 
'''il'^R = 4/9 is relatively large, hence the pacemakers are large in size and few 
in number (usually 1 — 3). In some cases, the scheme is able to remove all 
pacemakers after one application with vaccinating less than 10% of the popula- 
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tion. However, as the "pacemakers" compete with one another, there are cases 
where only the main pacemaker(s) is removed and secondary pacemaker(s) may 
appear in the next period. Eradication then requires a second application of 
the vaccination in the area of the remaining pacemaker(s). As shown in figure 
[l];,d, in such cases it typically requires vaccinating a total of some 18% of all 
individuals over both applications to bring the disease to total extinction. 

For the same model with syphilis parameters, the ratio ti/tr varies in a 
simulation within the range {1/8, 1/16}. As this ratio is small and variable, 
the pacemakers generated are more numerous (usually 2 — 5) and smaller in size. 
The additional model realism (20% of the population not gaining any immunity, 
birth/death process) makes the pacemakers "less circular" in shape (see Fig.|2^). 
As a consequence, it is necessary to vaccinate more areas, although each area 
is smaller in size. Nevertheless, the vaccination scheme generally eradicates the 
disease in a single application, requiring vaccination of 14% — 22% of individuals 
(see Fig. [2]) . As in practice it is difficult to obtain full coverage when vaccinating 
an entire population, or even a specific targeted group, to add realism (and lower 
vaccination rates) we vaccinated in Fig. [2] only an average of 86% of each of the 
identified areas, enclosing pacemakers. In more detail, the algorithm vaccinates 
up to 98% in the core of the pacemaker, where the high clustering of repeatedly 
infected individuals reside and as low as 60% (randomly chosen) in the outskirts 
of the pacemaker. 

It is of interest to examine the effects of simple random vaccination of the 
population. For the theoretical parameter values (a) based on [TP] and used 
in Fig. [l] the random scheme is only successful in eradicating the disease after 
vaccinating at least 80% of the population. For the syphilis parameter values, 
however, the random vaccination threshold is 43% of the population (this is 
somewhat similar to the results presented in [23] for an SIR model). Hence, a 
further gain can be achieved by combining random vaccination with the targeted 
scheme above. This advantageously reduces the vaccination threshold to some 
10% — 15% of the population for the syphilis parameters. Instead of vaccinating 
the entire area surrounding the pacemaker it suffices to randomly vaccinate only 
60% of the area normally targeted. 

In reality, this control scheme may be implemented by vaccination, quaran- 
tine or a targeted education plan, depending on the disease and on the means 
available for it's control (for information on efforts to eliminate syphilis in the 
US see [H] and CDC reports [23]). The main advantage of the vaccination 
methods proposed here is that they avoid the usual practice of vaccinating a 
large proportion of the population. In addition vaccination is confined solely to 
relatively small and specified areas (see figures [lj:,d and|2|D,c). In practice, it is 
always preferable to vaccinate as small a group as possible, as vaccination al- 
ways carries a risk. Hence, it is advantageous to target only the relevant groups, 
already at risk. In the case of syphilis, a vaccine [30] is still under development, 
but it is feared to be of relatively high risk - so if at all, vaccinating only carefully 
targeted population already at risk (e.g., in proximity to core groups of highly 
active individuals [3TJ|5]) will be desirable if and when such a vaccine is avail- 
able. Note, however, that while other works refer to tracing infected individuals 
or the most connected individuals for applying a targeted vaccination scheme 
(e.g., [5^131113]), in the scheme proposed here no contact tracing is required - 
the pacemaker waves stem from the SIRS dynamics and the small world struc- 
ture in a natural and intrinsic way. Hence, these areas are easily identified as 
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Figure 3: A bifurcation diagram of the proportion of infectives / for typical 
syphilis parameter values (as in Fig. 2) as a function of shortcuts p. For small p 
there is an endemic equilibrium. After a bifurcation point the dynamics exhibit 
a limit cycle with radius varying with p. The diagram plots the maximum and 
minimum values of / on the cycle. After a second bifurcation value of p the 
disease goes to extinction (due to synchronization). 



small areas where the infection appears repeatedly (see Figs. [T] and [2|. 

1.3 Synchronization and disease extinction 

Worthy of comment is the model's behavior for larger values of p, typically p > 
0.1, outside the "small world" regime, and corresponding to high levels of sexual 
promiscuity in the case of STD's. Counterintuitively, the epidemic consistently 
dies out abruptly due to the appearance of large-scale synchronized epidemics 
[ini UHl 123 - a well known cause of disease extinction. The synchronization 
manifests with the formation of large spatial aggregations of infected individuals. 
Upon recovery, these infectives gain temporary immunity for a lengthy time 
period. Thus the areas that once contained aggregations of infectives, become 
exhausted of susceptibles and there is no possibility for an epidemic to sustain 
- it soon dies out. 

Although it might at first seem unusual, the model implies that for society 
at large, grand sexual promiscuity has the potential to eliminate STD's such 
as syphilis altogether after about two decades of consistent behavior (see Fig. [s] 
and Fig. ^c)). The same would be true for populations in which individuals 
consistently have few proximate sexual partners. The most conducive conditions 
for the persistence of such STD's, appears to be the small world structure similar 
to the varying manifestation of sexual promiscuity seen in western society over 
the last centuries. 
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1.4 The impact of short-cuts 

The effect of the parameter p - the proportion of short-cuts, on the disease 
dynamics, may be assessed from the bifurcation diagram in Fig. |3] The figure 
plots the range in the number of infectives (maximum and minimum values) 
for any given p, when the model is run using the standard syphilis parameters. 
The following bifurcation scenario takes place: for p = the disease goes to 
extinction; for small p < 0.001 an endemic equilibrium is reached in which there 
is a relatively small proportion of infectives < /* <C 1; for 0.001 < p < 
0.13 (approximately) there is a limit cycle of radius depending on p and hence 
noticeable oscillations in the number of infectives; for p > 0.13 the disease goes 
to extinction due to a synchronization effect. The "limit cycle region" is the 
region where pacemaker centers develop and within it lies the region where the 
targeted vaccination scheme works very effectively. This region is an open strip 
of p values in the "small world" regime. Note that both for p = and for large 
p the disease rapidly goes to extinction. Thus, despite the period of temporary 
immunity built into this model, oscillations in I vanish for either very small or 
relatively large values of p (i.e., outside the small world regime). 

2 Difference equation model 

We formulate the following difference equation model to help gain insights into 
the network model's dynamics. Let St, It and Rt be the proportion of suscepti- 
ble, infective and recovered individuals in a large population at time t. Again, 
let Tj be the time period an individual remains infectious and tjj the period an 
individual remains immune. If we assume that r/ = 1, as the case for syphilis, 
the proportion of recovered individuals can be described by the sum: TiJ^^It-i, 
and thus: 

St = l-It^ S£T (2) 
where tq — tj + r^. The above model formulation is well known (see e.g., 
[35]), but is extended as follows. We suppose each individual has on average 
K connections including those to nearest neighbors and short-cuts. For the 
typical individual, denote by /C™^^ the number of nearest neighbors that are 
infected at time t and denote by Kl'^^^ the number of short-cut links that point 
to infected individuals. Then: 

= It{l-P)K, (3) 
= hpK. (4) 

Set Qp as the probability of being infected by a short-cut link and set qx as 
the probability of being infected by a nearest neighbor. In practice qp > qx, 
because of the important role short-cuts play in spreading the epidemic through 
the network. For example, simulations of the lattice model (Eqn. [T]) under 
syphilis parameters show that qp ~ Aqx- For the model parameterized with 
the theoretical values taken from |19j . qp ~ lOg/f . Incorporating this important 
observation in the difference equation leads to the following model of the SIRS 
dynamics on a complex contact network: 

/i+i = (1 - (1 - gK)^""(l - <7p)'^'"')5t, (5) 
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Figure 4: Time series of the difference equation model ^ for different p values 
close to bifurcation points; qK = 0.11, qp = 0.65, K = 12. 



where g/f < (7p < 1. Note that equation ^ captures both the time delay 
dynamics resulting from the temporary immunity of the disease and the social 
effects of the short-cuts. 

The dynamics of the model (|5| depend on p in a manner that is very similar 
to the more complex contact network model ([T|). Fig. |4] shows results based 
on setting qj^ = 0.11 and qp ~ 0.65. For p = 0, a small endemic equilibrium 
is reached (see FigJlFa)), for p in the small world range sustained oscillations 
arise (see Fig. 4(b)) |4p3)), and for large p the disease is eradicated (see Fig.Qc)) 
due to a synchronization effect. Comparing figures [2}i and 4(b) gb) - one can 
see exactly the same type of dynamics with similar period of about 10 years, as 
observed in syphilis datasets [231 E] • 

2.1 Stability analysis 

The equilibrium solutions of model (|5| are found by solving the following equa- 
tion for /* : 

r = (1 - (1 - (7A')'*(^-^)^'(l - 9p)''^^)(l - To/*). (6) 

It is easily seen that the infection-free equilibrium /* = is always a solution. A 
stability analysis (based on linearizion of Eqn. [5| reveals that the infection-free 
equilibrium is stable when the following inequality holds: 

Ro^~K[{l~p)lli{l-qK)+pHl~qp)]<l. (7) 

Hence, for the given model parameters and relevant K values (K = 8 or 
12), the infection free equilibrium is stable either only for extremely small p 
values {K = 8), or never {K = 12). This is visualized in Fig. [sjja) which 
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Figure 5: Bifurcation diagrams of the difference equation [s]; tq — 11, ~ 0.11. 
(a) = 0.65. Infection free equilibrium (/* = 0) is stable for parameters below 
solid curve along which i?o = 1- Dashed curve is the approximation based on 
estimating i?o (see text), (b) Hopf bifurcation curves. Endemic equilibrium is 
stable for parameters below indicated curves and a stable limit cycle is born 
when these curves are crossed (solid curve: K = 12, dashed curve: ii' = 8). 



plots the bifurcation curve of the infection free equilibrium as a function of the 
parameters K and p. The infection free equilibrium is stable for all values of 
parameters below the (lower solid) curve where i?o < 1 and unstable otherwise 
since i?o > 1. For if < 4, the infection free equilibrium is stable in a wide strip 
of the parameter plane, containing the small world regime (as visible in Fig. 5a). 

A similar result is obtained by estimating the reproductive number i?Oi ap- 
proximating it as the average number of secondary infectives produced by a 
typical infective individual in a sea of susceptibles. In the case here, where a 
single node may infect only those nodes it is linked to, the number of secondary 
infections may be approximated as: 

i?o«i^(<ZK(l-p)+'ZpP). (8) 

The above estimate gives a good approximation of the exact condition Q as 
seen in in Fig. [5]^ a) (dashed curve). 

For most parameter values, the difference equation ([s]) has a second endemic 
equilibrium in which /* > 0. The bifurcation curves describing it's stability in 
the (p, gp) parameter plane are shown in Fig. [sj^b) , for A' = 8 (upper curve) and 
K = Yl (lower curve). This curve indicates a Hopf bifurcation in the {p,qp) 
plane, where all other parameters are kept fixed. See Appendixl for technical 
details on its calculation. The endemic equilibrium exists for all parameter 
values for which the Hopf bifurcation curve exists. Below the bifurcation curve 
the endemic equilibrium is stable and is unstable above it, where a stable limit 
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cycle exists. The stable limit cycle thus coexists with the two unstable equilibria 
- infection free and endemic. 

Hence, the bifurcation scenario, for syphilis parameter values, for K — 8, is 
as follows. For < p < 0.0091, the infection free equilibrium is stable and the 
disease goes extinct; for 0.0091 < p ^ 0.079 the infection free equilibrium loses 
it's stability and an endemic equilibrium is born; for p > 0.079 the endemic 
equilibrium loses it's stability through a Hopf bifurcation and a stable limit 
cycle is born. See the upper dashed Hopf bifurcation curve in Fig.|5];b). In the 
case of K — 12, the infection free equilibrium is unstable for all values of p. 
However, the endemic equilibrium is stable for < p < 0.011. When p > 0.011, 
both equilibria are unstable and instead a stable limit cycle is observed. See the 
lower solid Hopf bifurcation curve in Fig. [5|b). 

As in the network model, the disease goes completely extinct for relatively 
large p values (for K = 12, pc = 0.43 and for K = 8, Pc = 0.87), despite the fact 
that the infection free equilibrium is unstable (see Fig. Qc)). The extinction 
should be attributed to the synchronization effect that takes place for these 
high p values. That is, the dynamics are such that a large proportion of the 
population become infected together and proceed on to move to the recovered 
class together. The synchronization requires the initiation of a strong epidemic, 
implying that Rq must be greater than unity, which explains why the infection 
free equilibrium is unstable in this regime. Nevertheless the disease becomes 
extinct due to the synchronization effect. 



2.2 Vaccination 

A random vaccination scheme may be incorporated into the difference equa- 
tion model by replacing St in Eqn. ^ with (1 — v)St. The parameter v is the 
proportion of susceptibles vaccinated per time unit. Denote by Ve the thresh- 
old proportion of vaccinated needed for disease extinction. A simple algebraic 
expression can be obtained for the extinction threshold by linearizing Eqn. ([s]) 
about the infection free equilibrium. Then, by using equation ([t]), disease ex- 
tinction is reached if: 

V>Ve = l-^. (9) 

For the parameter values of Figs. |4] and [S] the vaccination threshold is = 
0.5297. Numerical simulations corroborate the existence of this threshold. 

As the difference equation does not give any information regarding spatial 
patterns, it is impossible to apply the spatially oriented targeted vaccination 
scheme described above. However, targeted vaccination schemes may neverthe- 
less be explored by differentiating between vaccinating nearest neighbors and 
short-cuts. This can be achieved by replacing the term It in equations ^ and 
([5]) with the term {1 — VK)It for nearest neighbors and {1 — Vp)It for short-cuts. 
The results reveal that the extinction threshold is very sensitive to and low- 
ered dramatically by the term Vp for vaccinating short-cuts, while the term for 
vaccinating nearest neighbors, vk, has little influence. Although according to 
the model the infection free equilibrium is never stable for the syphilis disease 
in a small world type society, it appears that targeted vaccination effectively 
reduces the number of the actual contacts of the key individuals, thereby re- 
ducing i?o < 1 in the vaccinated population. A detailed study of the targeted 
vaccination is left for future work. 
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3 Conclusion and discussion 



Two models for studying the dynamics of diseases with temporary immunity 
in complex population networks have been proposed - a lattice model, which 
incorporates spatial information, and a difference equation model, which allows 
an analytic approach. The study focuses on the example of syphilis epidemics, 
which is a representative STD targeted to be eliminated in the US with little suc- 
cess so far. The network model reveals that diseases with temporary immunity 
on a small world contact network exhibit periodicity and waves of epidemics, 
some of which become pacemaker centers. It is shown that by eliminating pace- 
makers through vaccination, the disease goes to extinction within 1 — 2 periods, 
where only about 20% of the population requires vaccination. This is in con- 
trast to standard vaccination programs that set out to achieve herd immunity 
by vaccinating over 80% of the population. Moreover, by treating only those 
individuals in high risk pacemaker areas, it minimizes the application of the vac- 
cination with its possible risks to the larger population. The difference equation 
model allows further investigation of the Hopf bifurcation lying at the heart 
of the pacemaker phenomena. The models presented here were constructed in 
accordance to the US syphilis datasets (see [231 [TS] ^^nd [TT). The two models 
complement each other, allowing a more profound view of the dynamics. 

This work addresses the controversy as to whether syphilis epidemics recur 
approximately every ten years due to the temporary immunity it endows to 
infected individuals (IHj) or due to changing patterns in social behavior f |23l 
I16j). As shown here, both factors are crucial for recurrent syphilis epidemics. 
Thus, for example, oscillations cannot occur outside the small world regime even 
in the presence of strong temporary immunity. For zero or very small p values 
(corresponding to none, or a very few, short-cut links), an epidemics cannot 
develop. Moreover, the analysis performed on the difference equation model 
reveals that if all individuals in a small world type population network have 
only a small number of contacts, the infection free equilibrium is stable. In 
addition, on the other side of the scale, it is pointed out that for large enough 
p outside the small world region (corresponding to many short-cut links) a 
synchronization effect takes place, eradicating the epidemics due to exhaustion 
of susceptibles. In contrast, a society whose social behavior approximates a small 
world network with moderate heterogeneous levels of promiscuity would sustain 
the periodic recurrences of the syphilis epidemics approximately every ten years. 

Finally, we conjecture that complete disease extinction is nevertheless achiev- 
able by a targeted vaccination scheme similar to that presented here. The tar- 
geting and vaccination of key individuals, effectively reduces Rq to less than 
unity in the vaccinated population, thereby leading to disease extinction. 

A Appendixl: Technical details for the differ- 
ence equation model 

Here we present the technical details of the stability analysis performed for the 
difference equation model ([s]) . Consider the S'usceptible - Infectious - i?ecovered 
- Susceptible (SIRS) dynamics. Assume r/ is the amount of time units an 
individual spends in the Infectious class and that th represents the time units 
an individual later spends in the Recovered class. As for syphilis tj — 1 (where 
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a time unit represents 6 months), the proportion of recovered individuals can 
be described by the sum: 'E^^^It-i- Thus, denoting by St the proportion of 
susceptibles in the population at time t, we obtain: 

St = I- It- ^Ti'it-., (10) 

where tq = t/ + th. Set /( to be the proportion of infectives at time t, denote 
by K the number of connections an individual has on average and let p be the 
proportion of short-cut links an individual has among it's K connections. Then, 
denote by if"^ ^ the number of infectious nearest neighbors an individual node 
has at time t and by K™ ^ the number of infectious contacts via short-cut links 
an individual node has among its K connections at time t. Then: 



Kr'"" = it{i-p)K, (11) 

k';'"" = ItpK. (12) 

Set Qp as the probability of being infected via a short-cut link and set qk as 
the probability of being infected by a nearest neighbor, where qk < Qp < ^■ 
Hence, the SIRS dynamics on a network with nearest neighbors and short-cut 
connections can be described by: 

it+i = (1 - (1 - gK)^""(i - gp)^"')5*, (13) 



The dynamics of (13 1 are at equilibrium for the solutions, /*, of the equation 

r = (1 - (1 - qKY'^'-'^^'il - 9p)''^^)(l - to/*). (14) 

It is easily seen that the infection free equilibrium /* = is always a solution and 
that an endemic equilibrium /* > is a solution for most parameters relevant 
for syphilis. 

Substituting Jt — It ^ I* and linearizing Eqn.(13l about /*, we obtain: 
Jt+i = (l-<Z/f)^'-^^'^''(l-<?p)^^'* (i?o(l - To) Jt + S[l-i/t_,)-S[l-i/i_,+/i.o.t., 
where, 

i?0 - -K{{\ - p) ln(l - qK) + pln(l - qp)). 
Now substitute Jt = J{)\^ to obtain the characteristic polynomial, 

+ aA^°-i + I3X^°-^ + --- + PX + P^0, (15) 

where: 



a 



(1 - g/f)('-''^^''(l - qpr""'' (i?o(Tor - 1) - 1) + 1, (16) 
f3 = 1 - (1 - qK)^'-"^"'' {I - qpf"'' ■ 

The stability of the infection free equilibrium is derived from substituting /* = 
and requiring that |A| = Rq < 1. The stability analysis of the infection free 
equilibrium is presented in the main text (and Fig. [5ja)). Here we provide 
details regarding the endemic equilibrium stability and the calculation of the 
Hopf bifurcation curve(s), presented in Fig.jsjin the main text. This calculation 
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is inspired by a calculation of the Hopf bifurcation of a mean field difference 
equation model in |33j . 

Assume such a bifurcation exists and substitute A = e"^ (as stability changes 



at |A| = 1) into Eqn.(15l. Separating real and imaginary parts we obtain two 
equations: 

a — — cos((/)) — cot ^"^^ sin((^), (17) 
2csc(^)sin(|)sin(<?!)) 



/3 = 

COS 



\ ^ I \ / 

.s(^)-cos(^)- 



Substituting Eqns.(16) into Eqns.(17l and adding Eqn.(14l results in a system 
of three equations in the seven variables: qk, Qp^P, K, I* ,tq, (j). Some of these 
parameters can be fixed to values relevant for syphilis: To = 11, = 8 or 12 and 
g/f = 0.11. (f) can be viewed as the frequency of the periodic solution emerging 
at the bifurcation, where the period of the limit cycle (when exists) is T w 
As for the syphilis parameter values the period of oscillation is T w 20 time 
units, where 1 time unit = 0.5?; = 6 months, we fix (j) — 0.3. Now, with these 
fixed parameter values, we use the Newton method to solve Eqns.([T7| and ( [l4| , 
using the syphilis parameter values as the initial guess, once for K ~ 8 (dashed 
curve m Fig. [Sjb)) and once for K ~ 12 (solid curve in Fig.jsjb)). The results 
are plotted in the {p, Qp) plane (see Fig.jsf^b) in the main text). InFig.[5];b), the 
two bifurcation curves (one ioi K = 8 - dashed and one ioi K = 12 - solid) are 
presented, where below each curve an endemic equilibrium /* > is stable. At 
the bifurcation curve the equilibrium loses its stability and a stable limit cycle 
is born so that above the curve a stable limit cycle coexists with two repelling 
(unstable) equilibria - the infection free and an endemic. 
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